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ABSTRACT 


This  paper  reports  on  research  activities  developed  during  the  first  year  of  this  project.  The  main  tasks  accomplished 
and  reported  in  this  paper  are:  a)  the  development  of  new  computational  equation  of  state  (EOS)  for  granite.  One¬ 
dimensional  calculations  using  this  EOS  captures  the  first  principles  physics  of  the  near  source  region  and  identifies 
the  critical  material  properties  for  cavity  dynamics.  Cavity  dynamics  are  an  important  consideration  for  methods  that 
are  used  to  determine  accurate  estimates  of  yield  under  varying  emplacement  conditions;  and  b)  the  extension  of  the 
newly  developed  material  models  to  the  analysis  of  scaled  depth  of  burial  and  free  surface  effects  in  a  2D 
heterogeneous  structure.  Efforts  in  this  area  are  focused  on  a  better  understanding  of  the  source  stress-cage  region  as 
well  as  free  surface  Rayleigh  and  shear  wave  generation. 

A  strong  motion  code  for  uniform  source  region  structures  was  used  to  investigate  the  dependence  of  cavity 
dynamics  and  final  radius  on  material  properties  (Young’s  modulus,  shear  modulus,  gas  porosity,  overburden,  etc.). 
The  material  model  developed  was  obtained  by  taking  the  PILEDRIVER  and  the  HARDHAT  nuclear  test  events  as 
the  main  design  references.  The  following  aspects  of  the  problem  were  identified  as  the  driving  points  when 
developing  the  material  model:  velocity  profiles  at  given  stations  (near  field),  source  modeling  alternatives  (iron  pill 
vs.  ideal  gas  vs.  Hydses/SESAME),  energy  partition  after  the  shot,  peak  velocity  and  peak  displacement  attenuation 
profiles  and  cavity  size  as  a  function  of  the  depth  of  burial.  Previous  attempts  made  with  existing  material  models 
failed  to  comply  with  one  or  more  of  the  main  aspects  of  the  problem.  Because  of  this,  a  Tillotson  type  of  equation 
of  state  with  an  improved  handling  of  the  porosity  evolution  (crushing  and  bulking)  was  implemented,  along  with  a 
constant  yield  surface  strength  model. 

The  material  model  developed  is  used  in  a  set  of  2D,  axially  symmetric  simulations  aimed  at  studying  the  effects  of 
key  factors,  such  as:  the  variation  of  the  scaled  depth  of  burials  (SDoB)  and  the  effects  of  free  surface  topography  on 
the  cavity  dynamics  and  shear  wave  generation.  A  first  comparison  of  the  obtained  numerical  results  against  the 
empirical  prediction  equations  shows  that  the  evolution  of  the  calculated  final  radius  of  the  cavity  with  the  depth  of 
burial  is  relatively  closer  to  the  predictions  of  Heard’s  scaling  relationship  in  Mueller  and  Murphy  (1971)  than  to 
Denny  and  Johnson’s  (1991).  Further  research  needs  to  be  conducted  in  order  to  fully  understand  the  cause  of  the 
discrepancies  between  the  numerical  and  the  statistical  models.  Also,  the  cavity  growth,  the  predicted  spall  regions, 
and  the  seismic  radiation  associated  with  the  free  surface  effects  are  presented. 
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OBJECTIVES 

The  NNSA  Nuclear  Explosion  Monitoring  Research  and  Development  Roadmaps  (2010  these  Proceedings)  for 
Ground-Based  Systems  specifies  identification  of  the  sources  of  waveform  signals  as  an  important  research  and 
development  effort.  Source  identification  is  the  focus  of  this  project.  Understanding  the  physical  basis  of  seismic 
wave  generation  is  recognized  as  the  key  to  advancing  our  ability  to  understand  and  usefully  model  the  explosion 
source. 

The  first  main  objective  of  the  current  research  is  to  undertake  a  series  of  computer  simulations  to  study  the 
evolution  of  the  cavity  radius  under  different  conditions.  This  phase  of  computer  simulations  will  employ  LANL 
strong  motion  codes  (CASH)  for  uniform  source  region  structures  (i.e.,  ID  isotropic  velocity  models)  to  investigate 
the  dependence  of  cavity  dynamics  and  final  radius  on  material  properties,  such  as  Young’s  and  shear  moduli,  gas 
porosity,  overburden  (pgh),  regional  stresses,  etc.  The  calculations  were  done  using  a  material  having  some 
similarity  to  the  granodiorite  that  exists  at  the  Democratic  People’s  Republic  of  Korea  (DPRK)  test  site. 

The  second  main  objective  of  the  current  work  is  to  perform  2D  axially  symmetric  computer  simulations  that  will 
introduce  2D  sub-surface  heterogeneity  and  free  surface  topography  into  the  models.  Here  the  focus  will  be  to 
identify  the  impact  of  certain  idealized  source  region  structures  on  the  dynamics  of  cavity  growth  and  on  the 
radiated  seismic  energy.  These  2D  structures  will  be  introduced  into  models  of  the  source  region  in  controlled  ways 
in  order  to  isolate  the  effects  these  structures  have  on  the  close-in  phenomenology  and  to  identify  features  in  the 
cavity  growth  and  seismic  radiation  associated  with  those  effects. 

RESEARCH  ACCOMPLISHED 

Modeling  an  underground  nuclear  explosion  taking  into  account  the  propagation  of  the  wave  across  the  elastic 
region  comprises  a  wide  range  of  physical  and  thermodynamic  phenomena.  Four  different  zones  can  be  identified 
according  to  the  material  behavior.  The  first  zone  is  called  the  source  region,  where  the  rock  and  the  device  material 
are  vaporized  due  to  the  huge  amount  of  energy  released  from  the  explosion.  As  the  wave  created  by  the  explosion 
propagates  outwards  the  energy  dissipates  due  to  mechanical  and  geometric  factors.  At  a  certain  range  the  amount  of 
energy  is  no  longer  enough  to  vaporize  the  rock  but  it  is  sufficient  to  melt  it.  This  melted  region  defines  the  second 
zone  that  is  called  the  cavity  region.  As  the  wave  propagates  further  the  rock  material  does  not  melt  anymore,  but  it 
undergoes  plastic  deformations.  This  zone  is  called  the  plastic  zone.  The  commonly  known  shatter  zone  is  the  first 
part  (closer  to  the  source)  of  the  plastic  zone.  At  a  given  point  located  relatively  far  from  the  source  the  material 
starts  behaving  elastically,  i.e.  the  elastic  zone. 

In  order  to  properly  account  for  the  wave  propagation  across  all  the  zones  and  also  for  the  correct  energy  deposition 
in  the  elastic  region  of  the  rock,  a  series  of  key  constraints  concerning  each  one  of  the  zones  must  be  identified. 

These  constraints  are  obtained  from  the  available  experimental  data  and  are  to  be  employed  as  design  parameters  for 
the  purpose  of  the  development  of  the  appropriate  material  model.  In  this  work,  the  constraints  selected  are:  velocity 
and  displacement  profiles  at  given  stations,  energy  partition  after  the  shot  (for  each  one  of  the  zones  defined  above), 
peak  velocity  and  peak  displacement  attenuation  profiles  and  cavity  size  as  a  function  of  the  depth  of  burial. 

The  first  part  of  the  modeling  effort  comprised  the  development  of  the  material  models  to  be  used  in  the  different 
zones.  Regarding  the  source  region,  the  two  most  common  approaches  taken  in  the  past  are:  iron  gas  model  and  the 
bubble  model.  The  iron  gas  model  assumes  that  after  the  explosion  an  iron  gas  is  formed  occupying  a  relatively 
small  spherical  region  around  the  shot  point  (approximately  1.0  m  of  radius  for  a  10  kt  source,  [Schroeder,  1974]). 
The  initial  energy  of  the  explosion  is  therefore  uniformly  distributed  over  this  spherical  region.  The  bubble  model, 
on  the  other  hand,  considers  that  the  energy  of  the  shot  is  uniformly  distributed  over  a  spherical  region  with  radius 
equal  to  the  scaled  vaporization  radius,  r  =  2.0  FT1'3  [Schroeder,  1974],  where  W  is  the  yield  of  the  explosion, 
measured  in  kilotons  and  r  is  the  radius  of  the  spherical  region  in  meters.  Different  materials  have  been  used  in  the 
past  to  model  the  behavior  of  the  “bubble”  or  source  region.  The  two  most  common  materials  employed  for  this 
purpose  are:  ideal  gas  [Antoun  et  al.,  2001]  with  a  density  equal  to  the  density  of  the  rock  under  consideration,  and  a 
tabular  equation  of  state  (EOS),  i.e.,  Hydses  type  of  EOS. 

The  three  approaches  for  modeling  the  source  region  that  were  described  in  the  previous  paragraph  were  tested 
during  the  development  of  the  project.  However,  none  of  these  source  models  was  able  to  fulfill  simultaneously  all 
of  the  design  constraints  mentioned  before.  The  results  obtained  with  the  bubble  model  with  an  ideal  gas  provided  a 
very  good  match  with  the  experimental  velocity,  displacement,  peak  velocity  and  peak  displacement  attenuation 
profiles.  However,  it  over  predicted  the  final  size  of  the  cavity  and  also  the  energy  distribution  at  the  end  of  the 
simulation  was  not  realistic,  i.e.  the  amount  of  energy  remaining  in  the  source  region  was  only  4%  of  the  total 
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energy  of  the  explosion.  It  is  worth  noting  that,  depending  on  the  size  of  the  “bubble”,  the  amount  of  energy  needed 
to  vaporize  the  rock  contained  inside  it  ranges  from  27%  to  35%  of  the  total  energy  of  the  explosion.  On  the  other 
hand,  the  iron  pill  and  the  bubble  model  with  the  Hydses  granite  gave  reasonably  good  results  for  the  final  cavity 
size  and  for  the  energy  distribution  across  the  zones,  but  both  of  them  under  predicted  the  velocity,  displacement, 
peak  velocity  and  peak  displacement  attenuation  profiles. 

The  previously  described  source  model  approaches  arguably  do  not  fully  account  for  the  first  principle  physics  of 
the  shock  wave  continuum  from  the  cavity  expansion  all  the  way  up  to  the  final  energy  deposition  in  the  granite 
elastic  region,  i.e.  there  is  a  gap  between  the  near  source  region  and  the  near/far  field  region  dynamics.  The 
importance  for  future  studies  of  this  energy  deposition  and  the  proper  description  of  the  cavity  dynamics  combined 
with  the  shortcomings  of  the  previously  available  source  modeling  approaches  generated  the  necessity  of  finding  a 
different  approach  regarding  the  EOS  for  the  source  and  for  the  rest  of  the  granite.  Because  of  this,  the  Tillotson 
EOS  (Tillotson,  1962)  was  chosen  in  an  effort  to  find  a  bridge  for  the  existing  gap  between  the  near  source  physics 
phenomena  and  the  near/far  field  observation  points.  This  EOS  is  an  analytical  type  and  it  introduces  a  dependency 
on  the  specific  internal  energy  of  the  material,  which  makes  it  suitable  to  cover  all  the  zones  of  the  simulation  grid. 
The  phase  change  from  solid  to  vapor  is  also  taken  into  account.  According  to  the  state  of  the  material  (expansion  or 
compression)  and  to  the  level  of  specific  internal  energy  the  Tillotson  EOS  divides  the  p-v  space  into  four  different 
regions,  as  shown  in  Figure  1 .  The  Tillotson  EOS  utilizes  a  series  of  parameters  to  describe  the  behavior  of  the 
material  across  the  different  phases.  In  this  work,  the  parameter  set  used  was  based  on  Melosh,  (1989),  with  some 
modifications  that  were  introduced  in  order  to  adjust  the  material  model  to  the  lab  experiment  results.  The  source 
modeling  was  done  by  coupling  this  newly  implemented  EOS  with  the  bubble  model. 

The  material  model  is  completed  by  combining  the  EOS  of  the  rock  with  an  appropriate  strength  model.  The 
strength  model  describes  the  resistance  of  the  material  to  shear  deformation.  For  rock-like  materials,  the  shear 
strength  is  dependent  on  the  confinement  pressure.  The  main  component  of  the  strength  model  is  the  so  called  “yield 
surface”,  which  describes  how  the  shear  strength  changes  with  the  confinement  pressure.  The  yield  surface  used  in 
the  current  work  was  derived  from  the  one  proposed  by  Fossum  and  Brannon  (2004).  Some  modifications  were 
introduced  to  the  original  yield  surface  in  order  to  match  the  experimental  data  from  the  PILED  RIVER  and 
HARDHAT  events.  It  is  worth  noting  that  the  shear  strength  material  model  utilized  in  this  work  considers  that  the 
yield  surface  is  constant  during  the  whole  simulation,  i.e.  no  effects  of  strain  hardening  strain  softening,  high  energy 
states,  etc.  are  taken  into  account.  A  more  detailed  material  model  is  under  development,  where  all  of  these  effects 
are  to  be  taken  into  consideration.  The  yield  surface  used  in  the  current  work  is  shown  in  Figure  2. 

The  PILEDRIVER  and  HARDHAT  events  have  been  extensively  studied  in  the  past  and  very  good  matches  were 
obtained  between  the  numerical  simulation  results  and  the  experimental  readings  for  the  velocity,  displacement, 
peak  velocity  and  peak  displacement  profiles  (Antoun  et  al.,  2001;  Fossum  and  Brannon,  2004).  However,  to  the 
authors’  knowledge,  a  complete  analysis  of  these  events,  where  key  near  source  physics  parameters,  such  as  the  final 
cavity  size  and  final  energy  deposition  are  included  along  with  the  near/far  field  observational  data,  i.e.  peak 
attenuation  profiles  and  velocity  and  displacement  profiles  has  not  been  reported  yet.  Because  of  this,  the 
PILEDRIVER  and  HARDHAT  events  were  revisited  in  the  current  work. 

The  peak  displacement  and  peak  velocity  attenuation  with  the  scaled  range  obtained  from  the  numerical  simulations 
using  the  material  model  developed  are  compared  to  the  data  from  the  experiments  corresponding  to  the 
PILEDRIVER  and  HARDHAT  nuclear  events  (Rimer  et  ah,  1990)  in  Figure  3.  The  material  model  provides  a 
remarkable  good  match  for  the  attenuation  of  the  peak  displacements  of  the  HARDHAT  event.  In  the  case  of  the 
attenuation  of  the  peak  velocity,  the  numerical  results  over  estimate  the  experimental  data  for  HARDHAT.  This 
could  be  explained  by  the  fact  that  in  the  real  world  there  are  several  dissipation  factors  that  are  not  currently  being 
considered  in  the  material  model,  i.e.  weathered  material,  pre  existing  cracks  and  faults,  non  homogeneous  material 
formations,  etc.  An  interesting  feature  to  be  extracted  from  the  graphs  presented  in  Figure  3  is  the  clear  demarcation 
of  the  elastic  limit,  i.e.  the  point  where  the  material  ceases  to  behave  plastically  (yielding)  to  start  behaving 
completely  elastically.  For  the  material  model  developed  in  this  work,  this  threshold  is  located  at  a  scaled  distance  of 
approximately  200  m  /  kt1/3. 

The  PILEDRIVER  event  presents  a  trend  that  is  shifted  with  respect  to  the  HARDHAT  event,  even  though  both  of 
them  were  conducted  on  the  same  type  of  material  (NTS  granite).  Further  research  is  being  conducted  to  try  to 
explain  the  discrepancy  between  the  data  for  PILEDRIVER  and  HARDHAT. 

As  the  energy  from  the  source  propagates  across  the  computational  model  a  certain  amount  of  work  due  to 
compression  is  done  on  the  cells  of  the  model.  This  work,  which  is  essentially  p  dv,  where  p  is  the  pressure  and  dv  is 
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the  volume  change,  increases  the  specific  internal  energy  of  the  cell,  therefore  increasing  its  temperature.  When  the 
specific  internal  energy  reaches  a  given  value,  the  material  starts  to  melt.  As  it  was  already  pointed  out  before,  in 
this  work,  the  position  of  this  melting  point  was  selected  as  the  criterion  used  to  identify  where  the  boundary  of  the 
cavity  area  is  located,  i.e.,  the  last  melted  cell,  counted  by  starting  from  the  center  of  the  source,  is  located  and  the 
coordinate  of  its  center  is  taken  as  the  final  cavity  boundary  in  that  area. 

In  the  first  stage  of  this  project  the  evolution  of  the  final  cavity  size  with  the  depth  of  burial  was  investigated  with 
the  help  of  a  ID  spherical  computational  model.  A  constant  overburden  pressure,  corresponding  to  the  working 
point  depth,  was  applied  to  the  whole  model.  The  results  obtained  for  the  evolution  of  the  final  cavity  size  as  a 
function  of  the  depth  of  burial  are  shown  in  Figure  4. 

It  is  worth  noting  that  Heard’s  equation  for  the  cavity  radius  used  for  the  Mueller  and  Murphy  (1971)  curve  in 
Figure  4  was 

rc  =  16.3  W0'29  C  h~°A1  (1) 

where  W  is  the  yield  in  kilotons,  h  is  the  depth  of  burial  in  meters  and  C  is  a  material  dependent  constant  given  by 

C  =  h  p  ju  (2) 


where  E  and  p  are  the  Young’s  and  shear  moduli  respectively  and  p  is  the  material  density.  In  our  calculations 

C  =  1.4616  (3) 

The  equation  used  for  the  Denny  and  Johnson  [1991]  curve  in  Figure  4  was 


P 


147001T1/3 

0.3848  D0.2625 


(4) 


where  P  is  the  shear  wave  velocity  in  meters  per  second  and  P0  is  the  overburden  pressure  in  Pascal.  The  numerical 
values  for  the  material  properties  used  in  this  paper  are  given  by 

E  =  73.92  GPa 

p  =  2680  kg/m3  ^  ^ 

The  rest  of  the  variables  can  be  derived  from  these  two.  The  yield  was  kept  constant  at  one  kiloton. 

When  the  shot  is  done  at  deeply  over  buried  locations  the  calculated  cavity  size  is  very  close  to  the  predictions  of 
Mueller  and  Murphy  (1971).  As  the  depth  of  burial  gets  closer  to  the  nominal  value,  i.e.,  125  m  for  a  1  kiloton 
source,  the  calculated  final  cavity  size  gets  closer  to  the  predictions  done  by  Denny  and  Johnson  (1991). 

For  the  purpose  of  validation,  a  comparison  of  the  calculated  final  cavity  sizes  for  P1LEDRIVER  and  HARDHAT 
against  the  measured  quantities  is  presented  in  Table  1.  The  calculated  values  are  within  12%  of  the  measured. 

The  second  part  of  this  year’s  project  effort  involves  the  analysis  of  the  cavity  size  evolution  with  the  depth  of  burial 
utilizing  a  2D  axially  symmetric  computational  model.  These  2D  simulations  are  being  processed  at  the  time  this 
paper  is  being  prepared.  The  first  results  regarding  the  final  cavity  size  that  were  obtained  from  the  2D  simulations 
are  confirming  the  ID  values  shown  in  Figure  4. 

As  an  example  of  the  2D  results,  the  wave  propagation  pattern  0.3  s  after  the  shot  is  shown  in  Figure  5.  The  direct  P 
wave  resulting  from  the  explosion  is  shown,  followed  by  the  pP  wave,  which  is  a  result  of  the  bouncing  back  of  the 
direct  P  wave  against  the  free  surface.  The  shear  pS  wave  that  is  generated  as  a  result  of  the  interaction  between  the 
free  surface  and  the  main  wave  can  also  be  seen  in  Figure  5. 

The  computational  model  used  to  simulate  the  2D  examples  presented  in  this  work  does  not  include  the  possibility 
of  fracturing  or  spallation  inside  the  rock.  In  spite  of  this,  a  preliminary  prediction  of  the  spall  areas  can  still  be  made 
with  the  existing  computational  tool.  Each  material  has  a  certain  tensile  strength  that  defines  when  a  tensile  fracture 
would  start  to  develop.  Tensile  stresses  are  identified  by  negative  pressures  in  the  computational  model.  In  Figure  6 
four  different  snapshots  are  shown  where  the  spall  regions  are  marked.  In  this  case,  the  tensile  strength  of  the  rock 
was  assumed  to  be  equal  to  10  MPa. 
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Even  if  the  spall  areas  can  be  defined,  the  actual  influence  of  the  fracture  of  the  material  on  the  seismic  signals  is  not 
being  properly  captured,  i.e.  an  explicit  fracture  simulation  would  be  necessary  to  account  for  this  phenomenon. 
With  respect  to  this  matter  and  as  a  result  of  a  project  that  developed  over  the  last  years,  a  state  of  the  art  combined 
finite-discrete  element  methodology  (FEM-DEM)  (Munjiza,  2004)  was  recently  implemented  within  the  CASH 
hydrocode.  This  new  methodology  combines  the  ability  of  the  discrete  element  method  to  deal  with  a  large  number 
of  interacting  particles  while  keeping  the  finite  volume  ability  to  model  the  internal  deformation  of  the  material.  The 
areas  where  fracture  and/or  fragmentation  could  occur  are  modeled  as  a  collection  of  discrete  regions  or  elements 
connected  among  them  through  their  common  boundaries,  giving  a  continuum  like  behavior.  This  connection 
between  the  discrete  elements  is  described  through  the  appropriate  material  fracture  variables,  i.e.  fracture  surface 
energy,  maximum  normal  and  tangential  displacements,  etc.  In  this  way,  if  the  stress  state  across  the  discrete 
element  boundaries  is  such  that  the  fracture  toughness  of  the  material  is  not  exceeded,  the  material  behaves  as  a 
continuum.  However,  if  the  fracture  toughness  of  the  material  is  passed,  a  transition  from  continuum  to  discrete 
behavior  occurs,  with  a  series  of  fracture  patterns  as  a  result.  This  capability  will  be  very  helpful  in  the  future  of  the 
project  when  spallation  areas  are  to  be  defined,  and  their  effects  to  be  quantified. 

One  of  the  main  aspects  of  the  simulation  of  underground  explosion  is  to  be  able  to  properly  handle  extremely  high 
strain  rate  phenomena.  On  this  matter,  the  newly  implemented  FEM-DEM  methodology  was  tested  to  the  extreme, 
demonstrating  that  it  is  up  to  the  challenge  (Rougier  and  Knight,  2010a;  Rougier  and  Knight,  2010b;  Rougier  et  al., 
2010a). 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  source  model  approaches  used  in  the  past  arguably  do  not  fully  account  for  the  first  principle  physics  of  the 
shock  wave  continuum  from  the  cavity  expansion  all  the  way  up  to  the  final  energy  deposition  in  the  granite  elastic 
region,  i.e.  there  is  a  gap  between  the  near  source  region  and  the  near/far  field  region  dynamics.  In  an  effort  to  try  to 
close  this  gap  the  Tillotson  EOS  was  implemented.  The  verification  and  validation  of  this  EOS  was  done  by 
simulating  the  PILEDRIVER  and  HARD  HAT  events.  It  was  shown  that  the  material  model  developed  is  in  good 
agreement  with  the  HARD  HAT  experimental  data:  velocity  and  displacement  attenuation  profiles  as  well  as  final 
cavity  radius  prediction  and  expected  energy  distribution. 

This  material  model  was  used  to  predict  the  behavior  of  the  final  cavity  size  as  a  function  of  the  depth  of  burial  for  a 
constant  yield  of  one  kiloton.  The  results  obtained  showed  that  the  locus  of  the  calculated  final  cavity  size  is  located 
between  the  predicted  curves  of  Denny  and  Johnson  (1991)  and  Mueller  and  Murphy  (1971),  as  shown  in  Figure  4. 

It  can  be  seen  from  the  figure  that  Denny  and  Johnson  and  Heard’s  curves  from  Mueller  and  Murphy  have  a  stronger 
dependency  on  the  depth  of  burial  than  the  results  obtained  from  the  numerical  calculations. 

The  next  steps  on  the  project  will  be  concentrated  on  the  improvement  of  the  current  strength  model.  The  main 
modifications  that  are  planned  to  be  incorporated  are:  strain  softening,  strain  hardening  and  high  energy  effects.  A 
series  of  sensitivity  studies  are  to  be  conducted  after  the  improved  material  model  is  implemented.  These  studies 
will  mainly  concentrate  on  the  dependence  of  the  final  cavity  size  on  the  yield  of  the  explosion  and  on  the  main 
material  properties  (i.e.,  Poisson  ratio,  density,  “yield  surface”,  etc.)  These  studies  will  improve  our  understanding 
of  the  differences  between  the  statistical  and  calculated  hydrocode  results. 

The  capability  demonstrated  in  the  present  work  provides  a  valuable  tool  to  estimate  the  characteristics  of  unknown 
sources  that  have  been  tested  on  similar  materials,  like  the  case  of  the  DPRK  test  site. 

Further  studies  utilizing  the  above  defined  FEM-DEM  capabilities  coupled  with  the  newly  developed  equation  of 
state  will  enable  us  to  reach  a  new  plateau  of  understanding  of  how  the  seismic  shear  wave  generation  is  affected  by 
the  fracture  and  fragmentation  of  the  material  surrounding  the  source  area,  and  by  the  spallation  of  the  material 
located  above  the  working  point. 

It  is  worth  noting  that  this  new  EOS  and  material  model  approach  combined  with  the  first  principles  hydrocode  will 
be  applied  to  the  National  Center  for  Nuclear  Security  (NCNS)  2011  studies  at  the  Nevada  Test  Site  (Rougier  et  al. 
2010b). 
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Figure  1.  Tillotson  EOS.  Different  phases  in  the  p-v  plane.  -  Phase  I:  Compressed  states  for  all  values  of  e. 

Phase  II:  Expanded  state  with  e  <  eiv.  Phase  III:  Expanded  state  with  partial  vaporization  eiv  <e< 
ecv.  Phase  IV:  Expanded  state  with  complete  vaporization  ecv  <  e. 
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Pressure  [Pa] 


Figure  2.  Yield  surface  used  in  the  simulations  presented  in  the  current  work. 


Scaled  Range  (m/kt1/3) 


Scaled  Range  (m/kt1/3) 


Figure  3.  Material  model  validation:  Peak  radial  velocity  (left)  and  scaled  peak  radial  displacement  (right)  as 
a  function  of  the  scaled  range. 
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Evolution  of  the  Cavity  Radius  with  the  Scaled  Depth  of  Burial 


Figure  4.  Evolution  of  the  cavity  size  with  the  scaled  depth  of  burial.  Graphs  constructed  for  a  constant  yield 
of  1  kiloton. 
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Figure  6.  Prediction  of  the  spall  zone  for  a  lkt  charge  buried  at  nominal  DoB  (125  m). 


Table  1.  Comparison  of  measured  and  calculated  final  cavity  radius. 


PIFEDRIVER 

HARDHAT 

Measured  Final  Cavity  Radius  (m) 

44.5  19 

.4 

Calculated  Final  Cavity  Radius  (m) 

42.7  21 

.8 
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